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ABSTRACT 


The analvsis of stresses induced by contact between two bodies is inherently difficult 
because the size of the contact zone 1s unknown and constantly changing throughout 
loading. To overcome these difficulties, wo approximation methods have been devel- 
oped to determine the magnitude of contact stresses using the Ravleigh-Ritz method and 
the finite element method. Numerical optimization methods are employed to solve the 
contact problem. The solution techniques are compared to known analvtical solutions 
and shown to vield accurate results. An application of this approach to solving the 
contact problem is illustrated by examining the response of a clamped sandwich com- 
posite beam to low velocity impact. It was found that the maximum shear stress 1s in- 
sensitive to lamina thickness, however an increase in the contact laver thickness resulted 
in a reduction in interfacial shear stress. In addition, it was noted that a nonlinear 
bending stress distribution in the contact laver intensified as the thickness of this layer 
increased. This phenomenon was found to be localized to the region of contact. Finally, 
it was found that the compressive transverse normal stresses increased as the thickness 


of the contact lamina increased. 
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I. INTRODUCTION 


A. MOTIVATION 

Contact stresses occur when two bodies exert forces over limited contact regions. 
Examples of contact stress include meshing gear teeth, cam shaft and pushrod contact, 
rollers in plate forming operations, roller and ball bearings in contact with races, shaft 
and journal bearing contact, and plate-pin connections. 

Contact zones can be point, line, or surface in geometry. Because of the limited 
contact zone, the local stresses can be sufficiently high to be of major concern to the 
designer. Consequently, a thorough understanding of this phenomenon 1s essential. The 
first successful analytical solution to the contact problem between two spheres was 
solved in the late 19th century by Hertz. His solution can onlv be applied to simple 
geometries such as spheres, cvlinders, and flat plates. Because of these limitations, al- 
ternate solution techniques were needed to accommodate more complex geometries and 
boundary conditions. 

Unfortunately, the contact problem 1s difficult to study. The most significant diffi- 
culty is that the size of the contact zone is unknown and constantly changing throughout 
loading. Solutions are obtained by an iterative process. Consequently, the problem 1s 
highly nonlinear in its behavior. It is because of these difficulties that approximation 
and numerical methods are preferred in the solution of the contact problem. There are 
two general approaches used in the approximate techniques to solve the contact prob- 
lem. One class applies specific iterative procedures to solve a nonlinear system of 
equations that represent the contact state. For example, the contact condition can be 
simulated by the introduction of additional coupling terms into the system of equations. 
The other class constructs a functional that includes the contact body constraint. The 
functional is then minimized using specific numerical programming techniques. Al- 
though these two classes use different procedures, there are several methods to tincor- 
porate the contact condition into the problem formulation that are common to both 
classes. Examples of these include the penalty method and the augmented Lagrange 
multiplier method. These methods specify the manner in which the contact boundarv 
conditions are treated. 

The objective of this study has two parts. First, two approximate solution tech- 


niques will be developed to obtain the solution of the general contact problem. These 


techniques belong to the general class of methods that calculate a specific functional and 
then applies numerical programming methods to solve the contact problem. Second, 
these solution techniques will be verified by comparison with known analvtical solutions. 
With the solution techniques verified, the methods will be applied to study an actual’ 
contact problem. 

In order to accomplish the first objective, two models will be developed. The first 
method utlizes the Ravleigh-Ritz method to solve the contact problem. An assumed 
deformation field that satisfies the boundary conditions 1s carefully selected. Theory of 
elasticity relationships are applied to obtain an expression for the system's strain energy 
enabling calculation of the total potential. The minimization of the total potential en- 
ergy enables the calculation of the contact stresses at anv point in the body. The second 
model developed uses finite element analysis to accomplish the same objectives. The 
numerical minimization technique used in both cases is the augmented Lagrange multi- 
plier method. 

To accomplish the second objective, both methods will be verified by comparison 
with the Hertz solution of an infinitely long cvlinder in contact with a flat plane. With 
verification completed. the contact stresses resulting from low velocity impact between 


objects and composite sandwich materials will be studied. 


B. LITERATURE SURVEY 

As discussed previously, there are two general approaches to solving the contact 
problem. One approach uses special procedures to solve a nonlinear svstem of 
equations. The other creates a functional that includes the contact boundary con- 
straints. Although different, both approaches use similar mathematical methods to in- 
corporate the contact condition into the problem formulation. Two of these methods 
are the Lagrange method and the penalty method. A detailed explanation of the method 
of employing the Lagrange method in nonlinear finite element analysis was discussed bv 
Pian [Ref. |]. Alternately, the penalty finite element method was used by Cheng [Ref. 
2] to solve the multibody contact problem. The fundamental concept of the latter ap- 
proach 1s the transformation of a constrained problem into an unconstrained one. This 
is done by penalizing (i.e., increasing) the objective function for constraint violations. 
Although both methods are effective, there is substantial discussion on the limitations 
of both methods. Nour-Omid [Ref. 3] described the positive and negative aspects of 
both methods. The Lagrange method has been shown to be the more accurate method. 


However, its usage requires the introduction of additional unknowns thus increasing the 


tJ 


Il. FORMULATION OF THE CONTACT PROBLENI 


A. PRINCIPLE OF MINIMUM TOTAL POTENTIAL ENERGY 

As discussed in the introduction, the limited utility of the analytical solutions ne- 
cessitated the development of solution techniques capable of handling the nonlinear be- 
havior of the contact problem with complicated geometry and complex boundary 
conditions. This study intends to develop two numerical procedures to solve the contact 
problem. In short, the procedures will use different methods to obtain a functional, the 
system's total potential energy, and then use similar methods to obtain the equilibrium 
condition. Determination of the equilibrium position is made by application of the 
principle of minimum potential energy. With equilibrium established, contact stresses 
can be quantified. [n order to understand the details of this approach, the principle of 
virtual work and the principle of minimum potential energy must be discussed. 

Given a body in equilibrium, it 1s desired to describe the response of that body to 
infinitesimal displacements resulting from a system of forces. [f each particle in the body 
is described by some generalized coordinates, then the work resulting from these 
infinitesimal displacements is simply the product of the generalized forces acting on each 
particle and the particle’s displacement. However, if the particle is in equilibrium, the 
ioOikemist be Zero since the summation of forces in the x, y, and z directions 1s zero. 
The infinitesimal displacements and work in this example are referred to as virtual in 
nature. The fact that this work vanishes is referred to as the principle of virtual work. 

The virtual work discussed thus far can be subcategorized as virtual strain energy 
and virtual work done bv external forces. From the definition of strain energy, Virtual 
strain energy, OU, that results from virtual displacements can be calculated. Since this 
energy is viewed as energy against the bonds between elements, dU is a negative quan- 
titv. The work done by external forces is designated d+ and is simply the summation 
of the product of the external forces and the displacements of the generalized coordi- 
nates. 

Since the principle of virtual work states that the work done as a result of virtual 


displacements 1s Zero, 


é6W—dU=0 


Wi 


Alternately, this can be expressed as, 
= oC a0 


where [1 represents the system's total potential. 





The above equation illustrates the condition of minimum total potential of a svstem.. 
This is the foundation of the principle of minimum potenual energy. [Ref. 12: pp.. 
330-331] This principle states that in a condition of stable equilibrium, the svstem’s. 
total potential is stationary. Hence, determination of a system’s total potential energy 
and the minimization of that quantity will enable the calculation of the equilibrium po- 


sition. This is the basis for the numerical techniques developed in this study. 


B. CONTACT PROBLEM DESCRIPTION 

The objective of the solution techniques to be developed is to obtain a means of 
calculating the total potential energy and minimize it to determine the equilibrium con- 
dition. To accomplish this, two different models will be created to first solve a simple 
isotropic case, the solution of which is known. In this manner, our models can be vali- 
dated for use in more complicated arrangements. 

Consider contact between a cylinder and an infinite plane as shown in Figure lt. A 
cross section of the contact zone 1s shown in Figure 2. Let the width of the contact zone 
equal a distance of 2a. An analytical solution of this problem is available as a result of 
the work done by Hertz. In developing the solution methods, only a very limited region 
adjacent to the contact zone will be examined. The reason for this is that the contact 
phenomenon is a very local one. Figure 3 represents the analytical solution of the 
normal stresses resulting from this contact problem. [Ref. 13] As shown, the stresses 
in the foundation diminish very rapidly. The diagram shows that ao, becomes negligible 
at depths less than one half-contact zone (a) and a, diminishes significantly in less than 
3 half zones. Because of this, it is reasonable to assume that displacements beyond a 
very limited region are negligible in the strain energy calculation of the contact problem. 


A number of simplifying assumptons will be made for this study: 


|. As discussed above, displacements beyond a limited region are negligible in strain 
energy calculations. 


2. The foundation is an elastic isotropic material. The cylinder (roller) is rigid. 


system's degrees-of-freedom and computational time. The penalty method does not re- 
quire additional computational time, but it has been shown to result in less accurate 
solutions in satisfying the contact boundary conditions. Since the contact boundary 
conditions are exactly satisfied only when the penalty term goes to infinitv, the correct 
choice of the penalty parameter is the key to an acceptable solution. Guerra [Ref. 4] 
supported these claims. 


Because of the limitations of the Lagrange and penalty methods, Bishoff [Ref. 


CA 
=) 


advocated the use of the augmented Lagrange multiplier method to solve finite element 
contact problems. This method is favorable since it avoids the limitations of both 
methods discussed above. The augmented Lagrange multiplier method is similar to the 
penalty method in that the objective function is penalized for constraint violations. 
However, an additional multipller term is added so the optimum can be achieved bv a 
combination of both terms. As stated by Vanderplaats [Ref. 6: p. 141], the advantage 
of this method 1s that the penaltv term 1s not required to grow to infinitv to achieve exact 
constraint satisfaction. 

The augmented Lagrange multiplier method can be used in numerical programming 
techniques to find the optimum of anv functional. Pierre and Lowe [Ref. 7] provide a 
detailed analvsis of the programming techniques necessary in applving this method. 
Additionally, Vanderplaats [Ref. 6: pp. 140-147] provided an excellent discussion on the 
practical usage of this method in numerical techniques. Rothert et al. [Ref. $] used a 
numerical programming code based on this method to solve a nonlinear contact prob- 
lem. In this study, an existing numerical optimization routine uulizing the augmented 
Lagrange multiplier method will be used as an integral part of two solution methods 
developed to solve the contact problem. These methods will use different techniques to 
obtain the same functional. The augmented Lagrange multiplier method will then be 
used in a similar fashion to solve each problem. In each case, the optimization rouune 
will be used to determine a set of design variables that describes the contact state. 

Following the development and verification of the numerical procedures, this study 
will investigate the response of composite sandwich materials to low velocity impact. 
One of the common failure mode of low velocity impact is delamination. Joshi and Sun 
[Ref. 9] studied the impact response of a three laver cross-ply graphite epoxy laminate. 
A correlation was obtained between delamination cracks initiated experimentally and 
maximum shear stress points determined numerically. Sun and Rechak [Ref. 10] fol- 
lowed up these findings and found that the introduction of adhesive lavers between 


laminae reduced the shear stress distribution thus reducing delamination. Choi, Wang, 


and Chang (Ref. 11] studied the effects of laminae orientation, ply thickness, and 
stacking sequence on impact damage of graphite’epoxy composites. It was determined 
that stacking sequence affects impact damage more than laminae thickness variations. 
Much of the previous work has focused on the behavior of the graphite epoxy laminate. 
I{owever, there 1s currently interest in the development of turbine blades constructed of 
sandwich composites. [tis therefore beneficial to investigate the response of composite 


sandwich materials to low velocity impact. 





Figure 1. 


Figure 2. 


Roller-foundation assembly 


Contact zone cross section 


2a 


3. Deformations normal to the cross-section are negligible, hence a condition of plane’ 
strain eXists. 


4. The roller is subjected to a vertical distributed load. 


The roller-foundation contact js frictionless. 


Can 


An important restriction on the minimization problem will be that one body will be 
prohibited from penetrating into the other body. This mav seem like an obvious re- 
striction, however a method of mathematically stating this restriction must be discussed. 
Figure 4 shows the deformed and undeformed contact zone. Let 6 represent the de- 
flection of the roller due to an external force and v(.x,y) represent the vertical deflection 
of the foundation at any point (xy). At the point of Contact A, the condition of no 


intetference cam be EXptessed as, 


V(0,0) > 0 


At point B, this condition can be stated as, 
v(r sin 8,0) > 0 — r(l — cos 9) 


The latter condition can be specified at as manv points as necessarv to define this re- 


Striction. 


C. NUMERICAL OPTIMIZATION 
1. Optimization Fundamentals 

Before developing the models to be used in this study, one final area must be 
discussed. Once the total potential energy has been calculated, a method of minimizing 
it to find the equilibrium position must be employed. A study of the numerical opti- 
mization technique to be used 1s required. 

The technique being used in this study relies heavily on the methods of design 
Optimization. Design optimization is the utilization of mathematical techniques to 
minimize or maximize a particular value ‘to obtain an optimum solution. The method 
is ideally suited for design. A given design task may have an infinite number of sol- 
uuons. However, finding the best solution 1s a matter of the designer's experience and 
intuition. In the absence of significant experience in a particular field, finding this sol- 
ution may reduce to examining a range of possible solutions by trial and error. Opti- 
mization routines can be utilized to find this solution mathematically. 

Optimization problems can be constrained or unconstrained. For example, it 


may be desired to determine the minimum of the parabola, F(x) = (x — 5)? + 2. As seen 
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Figure 3. Analytical solutions for o, and o, 


in Figure 5, the minimum 1s clearly identified at point A. This is an example of an un- 


constrained problem. A constrained counterpart of this problem is: 


Minimize: F(x) 


Subject to: F(x) > .5x+4 


As seen in Figure 6, the minimum of the constrained problem 1s at point B. This is a 


simple illustration of constrained minimization. The contact problem is a constrained: 


minimization problem. 





Figure 4. Deformed contact zone 
The value to be minimized or maximized 1s referred to as the objective function. 
The parameters to be determined are referred to as design variables. The optimization 
process i$ an iterative one. The objective function is evaluated. The design vanables 
are varied thus obtaining a new objective function. If the difference 1s within a certain 
tolerance or meets a certain convergence criteria, the optimum has been obtained. If it 


is Out of tolerance, the cycle repeats. 


FIX) 





Figure 5. Unconstrained minimization 


F(X) 





Figure 6. Constrained minimization 


2. Augmented Lagrange Multiplier Method 
The contact problem belongs in the class of constrained problems. There are 
several techniques of solving constrained problems. The technique used in this study 
belongs to a class of solution techniques known as sequential unconstrained munimiza- 
tion techniques (SUMT). This class of techniques is designed for the general nonlinear 
problem. The fundamental concept behind this approach is that a constrained problem 
1s transformed into an unconstrained problem and the objective function 1s minimized 
using an unconstrained minimization technique. This transformation 1s accomplished 
by assessing a penalty to the objective function for constraint violations. For example, 
if the design variables are varied in such a way as to enter the region where a constraint 
is violated (1.e., the infeasible region), the objective function would be assessed a penalty 
feememencased jas hus, in Order to minimize the objective function, the desi@n variables 
remain within the region where no constraints are violated (1.e.. the feasible region). 
excl. (6: pp. 121-123] 
There are a number of methods within this class of techniques. They essentially 
differ in the wav in which penalties are assessed. The technique used in this studv 1s the 
augmented Lagrange multiplier (ALM) method. 


Given the constrained inequality optimization problem; 


Minimize: FLY) 


Subsect to. eUt) = 0,7 —"I2,....0 
The Augmented Lagrangian 1s defined as, 
mm 


A(X, Ap) = FU + > Ale) + 57) + ple) + 977} 


= 


~ 


where, 

X = vector containing the design variables 

/, = Lagrange multipliers 

p =penalty parameter 

s, =slack variables which convert inequality constraints to equality con- 
straints 
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The first two terms of -{(X,4.p) represent the Lagrangian. From the method of 
Lagrange multipliers, it is known that minimization of the Lagrangian represents opti- 
mality. Like simple problems where / 1s simply an additional unknown to obtain, in the 
ALM. method / is unknown. Hence, a mathematical routine based on constraint values 
is used to select and modify 4 for successive iterations. [Ref. 6: pp. 140-147] Initial 
selection of this term can have a significant impact of the problem’s convergence. 

‘\s discussed above, a penalty is assessed to the objective function for constraint 
violations. This feature is apparent by examining the last term. A constraint Violation 
results in a positive value for g(.\) thus resulting in an increase in 4. The value of p is 
a scaling term which is sequentially increased throughout optimization. This ensures 
that there is a balance between convergence and numerical conditioning. If p remains 
small, convergence may occur with major constraint violations. If p remains large, 
constraints will be satisfied at the expense of an ill-conditioned problem [Ref. 7: p. 
169]. As with the Lagrange term, selection of this term has a significant impact on the 
outcome of the problem. 

3. Optimizer and One-Dimensional Search 

Thus far, a procedure has been defined which has transformed a constrained 
minimization problem into an unconstrained one. This level of the optimization process 
is referred to as the optimization strategy. The formulation of the modified objective 
function via the augmented Lagrange multiplier method represents a kev portion of the 
Optimization process. However, there are additional parts of this process that require 
comment. 

With the modified objective function and a procedure for assessing penalties in 
place, a procedure for minimizing the objective function must be defined. This portion 
of the process is carried out by the ‘optimizer.’ The optimizer’s function is to system- 
atically alter the design variables in a manner that reduces the objective function rapidly. 
If .. represents a vector containing the design variables, the following process 1s used by 


the optimizer to alter X 
ae ae X(j-1) == KS; 
where, 
i= current iteration number 


S = search vector 


k = scalar representing distance traveled in direction S 
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In general, two processes must be accomplished to find the optimum. First, the search 
direction S must be determined by a systematic process. An example of this phase is the 
steepest descent method where the direction of steepest gradient is chosen. Second, the 
scalar kK must be determined such that the objective function is minimized as much as 
possible in the search direction of the current iteration. The latter phase is referred to 
as one-dimensional search. [Ref. 6: pp. 10-12] 

The optimizer used in this study is the variable metric method. Due to the 
complexity of this method. a discussion of their formulation 1s omitted. A detailed ac- 
count 1s available in Vanderplaats [Ref. 6: pp. 92-93]. The one-dimensional search 
routine used in this study is the golden section method with polynomial interpolation. 
The one-dimensional search portion of the process merely represents a systematic and 
efficient method of finding the minimum in the chosen search direction. A detailed ac- 
count of this approach is again available in Vanderplaats [Ref. 6: pp. 26-49]. 

4+. Convergence 

The final point to be discussed relevant to optimization fundamentals is that of 
convergence. Convergence criteria are utilized to idenufy the optimum solution and 
terminate calculations. There are a number of convergence criteria that can be uulized. 
The most obvious 1s absolute convergence where the objective functions from two suc- 
cessive iterations are compared. [f the difference between the values is within some 
prescribed limit, optinuzation is terminated. A second method signifving optimality for 
unconstrained problems is calculation of the gradient with respect to the design variable 
vector VY. [fthis value is approximately zero, optimality has been achieved. This method 
is called the Kuhn-Tucker conditions for unconstrained minimization. Kuhn- Tucker 
conditions are more involved for constrained problems. [Ref. 6: pp. 100-101] 

The Automated Design Synthesis (ADS) System used in this study uses both 
these termination criteria as Well as relative convergence. Relative convergence 1s similar 
to absolute convergence except normalized versions of the difference between successive 
iterations is calculated. Again, if a specified tolerance 1s achieved, optimization 1s ter- 


minated. 


Hl. APPROAIMATE SOLUTION TECHNIQUES 


A. RAYLEIGH-RITZ APPROACH 
1. Background 





The Ravleigh-Ritz method 1s a method of utilizing the theory of minimum po- 
tential energy to solve a given problem. The fundamental concept behind this method. 
is that a trial function that represents the deformation field is chosen in terms of un-: 
known constants. Second, the system's total potential energy is calculated in terms of, 
the trial function. Since the total potential is a minimum at equilibrium, minimization: 
enables determination of the unknown constants. The total potential is minimized by 
differentiating with respect to the unknown constants and equating to zero. The result 
is ‘n’ equations and ‘n’ unknowns, the trial function constants. [Ref. 12: pp. 335-336 ] 


The only requirement of the Ravleigh-Ritz method is that the trial function is 


kinematically admissible. A kinematically admissible solution 1s one that satisfies the . 
geometric boundary conditions of the svstem (1.e., deflection and slope). Other require- 
ments need not be satisfied. As an example, consider a simplv supported beam of length 
L with the origin at the left end of the beam. A kinematically admissible solution to 
describe the beam’s one dimensional deformation from vertical loading in the y direction 
1S, 


(x) = a, sin( =) 


ip 


where a, represents the coefficient to be determined. 


Deflection boundary conditions at x =0 and x = L have been satisfied. Obviously, an 
increased number of terms in the trial function will vield a far more accurate solution. 
A Fourier sine series would be a reasonable selection in this case. 

Although the Rayleigh-Ritz method does not stipulate numerous requirements 
on the trial function, sensible choices of trial functions will increase solution accuracy 
significantly. For example, consider the same simply supported beam with the origin 


now at the center. A kinematically admissible function 1s, 





Naturally however, due to the placement of the origin in this problem, an even function 
is a much for sensible selection for a trial function. A more appropriate selection would 
be, 


TX 


iE 


j(x) =a, cos 


The latter point concerning sensible selection of the trial function will be discussed in 
detail throughout this study. 
2. Application of the Rayleigh-Ritz Method to the Contact Problem. 

A brief overview of the method to be developed is in order. Application of the 
Ravleigh-Ritz method necessitates the selection of the appropriate trial function in terms 
of unknown coefficients. A discussion of the physical nature of this problem as well as 
the desired features of the trial solution 1s required. 

The theorv of elasticity relationships will be applied using the trial function to 
obtain the svstem’s strain energy in terms of unknown coefficients. The svstem’s total 
potential energy will then be minimized uulizing the optimization techniques discussed 
in Chapter II. The design variables are the unknown trial function coefficients. With 
the coefficients determined, the displacement is known for all points enabling the stress 
MoEGeseteniimeotapouecnout the body. Figure 7 is a flow chart of the procedure to be 
utilized. The post-processing procedure shown is simply the calculation of the stresses 
using the now determined coefficients. 

3. Trial Function Selection 

In order to choose an appropniate tnal function, it 1s necessary to have an 
understanding of the physical phenomena to be modeled. Consider the roller-foundation 
svstem shown in Figure |. Two trial functions are needed to model horizontal and ver- 
tical deformation. There are a number of important characteristics that should be in- 
herent within the trial functions. These are outlined below: 


|. As seen in Figure 3, o, and o, are equal and compressive at the point of contact. 
Additionally, o, decays more rapidly than oa, as distance from the contact point in- 
creases. 


th 


Taking the origin at the point of contact as shown in Figure 2, the u-deformation 
takes the form of an odd function. 


The v-deformation takes the form of an even function (i.e. symmetric about the 
origin and greatest at the origin). 


loa 


4. To satisfy the Rayleigh-Ritz method requirements, the boundary conditions must 
be satisfied. In this case, this requires u and v deformation to be Zero at the far 
boundaries. 
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Figure 7. Rayleigh-Ritz method applied to contact problem 


5. The characteristics outlined in item 4 satisfy the requirements of this method. 
However, since this study will be a stress analysis, it 1s also desired that the stresses 
also reflect the physical phenomena. Since o, and o, are functions of ¢, and e,. €, 
and ¢, should also exhibit certain characteristics. Referring to Figure 8 for dimen- 
sions and the coordinate svstem, it 1s desired that €, and c, equal zero atx}=L 2 and 
[it nis “waillsensure tiatestresses are zero at the boundaries. To satisfy the re- 
quirements of item | above, ¢€, and €, should be equal and negative at the point of 
contact and decrease in magnitude as the distance from the point of contact 1n- 
creases. 





Figure 8. Contact zone 


With the above guidelines in mind, the tral function can be selected. The trial 


functions chosen for this study are composed of a sernes of terms of the general form: 


u(xy) = a,(H — y)"( —x) (3.1) 
C L d . 4 
v(xy) = b(H —¥)"(S -2) (3.2) 


These expressions were carefully chosen and represent a compromise due to the difh- 
culties of satisfying all boundary conditions with the physical phenomena of this prob- 


lem. 
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From the above expressions, it is immediately obvious that the geometric 
boundary conditions at x=L 2 and y=H have been satisfied. This satisfies the re- 


quirements of the Rayleigh-Ritz method. In addition, there are a number of important 





characteristics that illustrate the advantage of this selection: 


l. €, and é€, are negative thus simulating a compressive environment in the vicinity of 
the point of contact. The importance of this is obvious. If normal strains were not 
negative, the resulting requirement would be for the coefficient to be less than zero 
to simulate compression. It is obvious that this would result in deformations op-: 
posite to that which was desired by the choice of the trial function. 


This selection for deformation fields has the important characteristic of decreasing: 
deformation as we move away from the point of contact. Also note that defor-: 
mation 1§ maximum at the point of contact. 


tv 


The exponents a.b.c, and d can be varied to simulate subsurface stress fields. Fom 
example, if the analvtical solution indicates a large v gradient for ao, at x=0, the’ 
objective would be to increase the rate at which e, decreases as the distance from 
the point of contact increases. This would be easily simulated by raising the value: 
ofa. If this change had a detrimental effect on the behavior of o,, the exponents: 
of the vertical deformation could be varied to restore the solution. 


or 


[t 1s important to note that this selection 1s not without compromise. The most’ 
significant limitation of the trial functions is with regard to the horizontal deformation 
u. Phvsically, it 1s expected the u(x) behave as an odd function as discussed above; 
However, in this selection of trial function, a positive value of u(x’) exists at the origin. 
This 1s contrary to the physical behavior of the problem and will lead to some error. 
However, considering that the magnitude of this deformation in the elastic range 1s 
small, this error 1s believed to be limited. Another consequence of this compromise 1s 
the existence of non-zero shear strain at x = 0. 

Another less severe limitation is a restriction on the order of the exponents in 
order to maintain zero stress at the boundaries. Since o, and ao, are combinations of é, 
and e,, both normal strains must be zero at the boundaries to ensure that stresses are 
zero at these locations. Since the normal strains are first derivatives, this requires that 
b and c are at least equal to 2. 

It is Worthwhile to note that most of the considerations discussed above far ex- 
ceed the requirements stipulated by the Rayleigh-Ritz method. The objective has been 
to utilize trial functions that closely match the physical nature of the problem in an ef- 
fort to maximize accuracy. 

The specific trial functions used in this study are listed below. The horizontal 


deformation was assumed to be; 


uxy) =) a,(H — yt & — x) (3.3) 


ies 


— 


where // and L represent the height and length of the bearing foundation, respectively. 


Summation was done for n equal 1 and 4. The vertical deformation was assumed to be; 


vixy) = > bal H —y)*( ee Se ier (3.4) 
i= . 


As discussed previously, manipulation of the exponents enables the trial function results 
to be matched with the analytical solution. As will be illustrated in the results, the ex- 
ponents chosen in the above functions achieve this goal sufficiently. 

With the deformations chosen, the stresses and strains can be calculated for use 


in the calculation of the foundation strain energy. These values are shown below: 


f=) —2a,(H— yt x) 3.5) 

a! 7 

th 2b (H < Ltn) . 

fy = = al — y)( 5 x) (3.6) 

i=1 

oo oa [C1 — ve, +ve,] (3.7a) 
* (1+ v)(1 — 2v) eee 
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where E and v are Young’s modulus and Poisson’s ratio, respectively. For the shear 


Stress and strain, 





oH = Y —(4 4+ n)a,(H — yl S = x) (3.8) 
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Cu CV 
La en arr (3.10) 
xy — Gy (3.11)) 


where G 1s the shear modulus. 


Lsing the above quantines, the strain energy U can be calculated. From the definition. 


- 


of strain energy applied to a two dimension deformation field, 


Le eres 
: 2 A 
C= | ; ey 1G Ete ey tata iy te (3a 
= aw 
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Because of svmmetry about the origin. strain energy can be calculated for half of the 
domain and doubled. With strain energy calculated, the total potential for the system 


can be foulmd trem. 


where, 


F = external force per unit length applied to the roller 


s 


0 


vertical distance traveled bpathe roller. 


The quantity Fo represents the work done by the roller on the bearing foundation. 


B. FINITE ELEMENT APPROACH 
1. Total Potential Derivation 

The finite element method can be employed to solve the contact problem. A 
finite element mesh can be constructed to approximate the behavior of an elastic foun- 
dation subjected to line contact loading from an nigid roller. The resultant interaction 
between the foundation and roller enables the calculation of the foundation's strain en- 
ergy and the subsequent calculation of the total potential energy. By again utlizing the 
optimization techniques discussed in Chapter II, the equilibrium position can be deter- 
mined. Thus the contact stresses can be calculated throughout the body. 

The objective is to derive a means of calculating the total potential of the system 


by application of the finite element technique. Total potential energy is defined as, 


(3.14) 


[l= U — Fo 


where, 
C = internal strain energy 


j= elem eonce per unit lenetiapplied to the roller 


0 = vertical distance traveled by the roller. 


The strain energy of the svstem can be found from, 


C= | a Omer + Oye, =e ep dee (Se 13) 
Q 
where {2 represents the problem domain. This can be expressed in matrix form as, 
ve | Liao fenta) 
ae 
where, 
T ; 
{e}° = {ey ey V xy! 
fe 
{a} a oe Oy ae 
Winethe element level; 
Uy | pak 
v=) Se ace (3.17) 
Q' 
Where 1 represents analysis of the :” element. 
The stress matrix can be expressed as, 
(3.18) 


{o} = [D]{e} 


where [D] represents the material property matrix. 


For a condition of plane strain, the stiffness matrix can be expressed as, 





l-\ v 0) 
pone é we. (3.19) 
(L+vj(1 — 2v) (lee | 
0 0 = 
for a plane stress condition, 
ee 0 
)| = ——_ (3.196) 
(1 0M ) 0) () | ~ v 
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The development of this technique will use linear triangular elements. The method, 


however, can be applied to any type element. For linear triangular element. the defor-- 


mations take the following form: 
u= Ayu + Hyu, + Ayu, (3.20% 
v= Hv, + Hyv, + Av; (3.219 
Where the shape functions //, are defined as: 
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A, = (ah — a) Oa 
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Hy = [xy — <3) © 3 =) + = 3) pars | 


| ; 
Ay = [Onn — 221) + O1 2) + 02 — ISG (3.24) 





where, 
x,y, = coordinates for node 


A = element area [Ref. 14 }. 


The strain matrix can be expressed as, 





= 0 
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Substituting equations (3.20) and (3.21) into equation (3.25) vields, 


CH mid CH 1d 
l 0 Cily 0 3 0 




















eis Cx Ox Vy 
CH, Gilde CH, ls 
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In abbreviated notation. equation (3.26) 1s expressed as, 
{e} = [B]{d} 
For the linear triangular element, [8B] reduces to, 
5) a a ii!) 22 
(B= | 0 y-m 0 m-y 0 
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Returning to the element strain energy calculation, equation (3.17), 


and substituting (3.18) into (3.17), 


c= | > fe} [D]{ehaV. 
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Substttuting equation (3.27) yields, 


u=| + (CB]{d}) TD Bi ddV 
a 


vata tayo | ar 
, 


t= = (TB) [DIB] And} (3.30) 


where ¢ = unit depth. Defining the element stiffness matrix [A] , 
[A] =(B](D][B) 4 (3.31) 


then the strain energy per element equals, 
: l aa 
C= > (a3[Alia} (3.39 


The element stiffness matrix can be expanded into the global stiffness matrix. 
With strain energy now determined, total potential can be determined from equation 
(3.14). 

2. Optimization and Static Condensation 

As with the previously developed model, the augmented Lagrange multiplier: 
method will be uulized to determine the equilibrium condition via the theorem of mini- 
mum potenual. The objective function is again the total potential. In this case, how- 
ever, the design variables are the nodal deformations, u, and v, for non-fixed nodes. 
Constraint equations are developed in a similar manner as discussed in Chapter II to: 
ensure that one body does not violate space occupied by the other. 

Since the nodal deformations are represented as the optimization design vari-. 
ables, the number of design variables will equal twice the number of non-fixed nodes. 
For a simple mesh, a direct application of this procedure will likely yield accurate results. 
However, it is known that the accuracy of the optimization routine declines as the: 
number of design variables increases. Hence, for complicated finite element meshes, 
solution accuracy will be adversely affected by the large number of design variables. 


Therefore a procedure must be adopted to eliminate the need for assigning design vari- 





ables to nodes where information 1s not necessary for evaluating a solution. This pro- 
cedure 1s Known as static condensation. 

Static condensation has been utilized by References 2? and 3 in an effort to re- 
duce computer computational time. The idea behind static condensation 1s the reor- 


ganization of the global stiffness matrix. A finite element problem can be expressed as, 
CK, = (F) (3.33) 


where, 
[A,] = the stiffness matrix 


{uw} = the deformation vector 


{F} = the force vector 


It is desired to reorganize this svstem of equations into the following: 


Ary Aya | fy 7 fF 
A>, Ax | U4 0 
The vector uw, contains the essential nodes while vector uw, contains non-essential nodes. 


Essential nodes are those Where boundary conditions are applied and nodes that are as- 


signed optimization design vanables. By matrix manipulation, 
“ —| P \y 32 
{uy} =— LA.) [Ag Jiu}. [3.55)) 


Therefore, the displacement vector can be expressed as, 


a - 4) = om ian) (3.36) 
— [A] (Kp, Jiu} — [Ky] [Aa] } | 


where I represents the identity matrix. 


Substituting this equation into the global counterpart of equation (3.32), 
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By defining the reduced stiffness matrix [K], the above reduces to, 


U => (uy) TAH} (3.38) 


As discussed at the beginning of this section, the ADS design variables are as-: 






signed as the honzontal and vertical deformations at all non-fixed nodes. With the in-- 
tegration of static condensation, design variable assignments are further restricted to) 
non-fixed, non-condensed nodes. The procedure 1s now in place for calculation of strain; 
energy and total potential energy. In summary, a flow chart of the solution procedure: 


utilized in this chapter is shown in Figure 9. 


Initialize Program | 


Define Optimization Design Variables 
Calculate Element Stiffness Matrix [k} 
Assemble Global Stiffness Matrix [{KG] 


Reduce Degree of Freedom by Static Condensation {KG} 


Assign Design Variables to {d)} 


Calculate Strain Energy U 






ue(a) *[KG){4) 
Alter Design 
Variables 


Calculate Total Potential 
Define Constraints 


Call Optimizer to Minimize 
Total Potential 


No 









Solution 
Converged? 
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Figure 9. Finite element method applied to contact problem 
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IV. RESULTS AND DISCUSSION 


A. PROCEDURE VALIDATIONS 
l. Ravleigh-Ritz Method Results 
[In Chapter III an approximation technique was developed to solve the contact! 
problem via the Rayleigh-Ritz method. As discussed, two trial functions that approxi-- 
mated the deformation field were selected in terms of unknown coefficients. The hori-- 
Zontal deformation was assumed to be; 
m 


u(x) = > a(H ay) ha i = x) 


n=! 


The vertical deformation was assumed to be: 


Blagg = » b,(H or a _ ae 
= = 


Lsing the above deformation fields, theory of elasticity relationships and the definition 
of strain energy were emploved to obtain an expression for the total potential energy of 
the svstem shown in Figure |. Numerical minimization techniques were then employed 
to determine the equilibrium condition and the contact stresses. 

To illustrate the application of this method, an isotropic material with the fol- 


lowing properties was selected: 


E = 200 GPa 
v=0.3 
G = 76.9 GPa 


As stated in Chapter II, this problem was selected for development of this technique 
because an analytical solution is available as a result of the work done be Hertz. It is 
desired to use this analytical solution to choose a roller size, load, and problem domain 
that can be used to accuratelv simulates the contact phenomenon. 

Contact stresses as well as the size of the surrounding region of influence are 


strongly affected by the size of the contact zone. Naturally, as the size of the contact 
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zone increases, the load 1s distributed over a larger area and contact stresses decrease. 
The extent of the affected subsurface zone also decreases. From the analytical solution, 
the contact Zone size is defined by the externally applied force, the material properties 
and the diameter of the roller [ Ref. 13 ]. Hence for a given material. the roller size and 
the external load define the contact zone size and the resulting stresses. 

Figure 2 defined the width of the contact zone as 2a. Figure 3 shows the ana- 
Ivtical solution of the contact problem. This figure shows the decrease of the normal 
stresses as a function of half-contact zones (a) away from the contact point. As shown, 
o, decreases more gradually than o,. Therefore the decay of oa, is the limiting factor in 
defining the domain bevond which strain energy contributions are negligible. From 
Figure 3, 1t 1s estimated that the contact phenomenon can be accurately modeled by 
examuning a region equal to approximately five half-zones (Sa). 

Since a numerical integration technique was used to perform the double inte- 
gration required by Equation 3.12, the dimensions Were selected for numerical conven- 
mence. Referring to Figure 8, height H and length L were selected as 1 and 2 meters, 
respectively. Due to the problem’s symmetry, half the foundation was analyzed. This 
enabled the double integration to be conducted between the limits of 0 and 1. Since the 
foundation height H has been set to 1 m, it 1s desired to have this distance equal to 5 
contact zones (5a) as described above. 

Lsing the analytical solution, a load and roller diameter were selected that cre- 
ated a contact zone such that the foundation height H was equal to a distance of Sa. 
The only additional restriction was that the resulting contact stresses remained within 
the elastic range of the material. Yield stress was assumed to be 300 \f{Pa. The load 


and roller radius combination used in this study are; 


Load: 90 MN 
‘Radius: 75 m 


The values were obtained using the analytical solutions found in Reference 13. The 
latter dimension is unrealistic. However, as stated above, it 1s a result of the selection 
of base dimensions in the interest of numerical convenience. This model simply repres- 
ents a scaled-up examination of the base material in a small region adjacent to the con- 
tact zone. 

Comparisons of the approximated contact stresses with the analytical solution 


along the axis of symmetry are shown in Figures 10 and 11. Both figures are normalized 
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graphs of the stress distribution. As shown in Figure 3, the maximum stress occurs at 
the point of contact. As shown, this stress is equal to o, and a, at the point of contact. 
Figure 10 1s a comparison for a, . This figure shows a stress distribution that closely 
approximates the analytical soluuon. Figure 11 shows the analytical and approximate 
stress distribution for o,. It is apparent that this method over approximates this stress. 
As stated in Chapter III, one of the benefits of this trial function is the ability to change 
the exponents to match analytical solutions with approximate solutions. <A brief expla- 
nauon of the choice of exponents used in this solution technique follows. 

While selecting exponents, it must be understood that o, and a, are both func- 
tions of horizontal and vertical deformations. Therefore changing the exponent of one! 
deformation to alter the stress in one direction will influence the behavior of the other. 
In the case of Figure I], it appears that a modification is required. An increase in the 
exponent of the v-portion of the horizontal deformation function seems appropriate to. 
increase the rate at which a, decavs. However, this action will have the undesirable ef- 
fect of decreasing the rate of decay of o,. It is possible to counter by decreasing the: 
exponents of the vertical deformation. However. as stated in Chapter III, in order to 
maintain a zero stress boundary condition at x=L 2 and y=H the exponents of all 
terms must be greater than or equal to 2. Thus a compromise must be reached) It @ 
believed that since a, is the dominant term, priority should be placed on ensuring a, is 
as close as possible to the analytical solution. Accordingly, a decision was made to ac-: 
cept the over estimation of o, as shown in Figure 11. In this case, the estimation of o,. 
IS COnSenatinie. 

The contours of the Rayleigh-Ritz solution are shown in Figures 12 through 15.. 
It has been determined that this method can accurately predict contact stresses resulting: 
from line contact between a roller and flat plane. Table | compares the approximated! 


contact stresses and the analytical results at the point of contact. 


Table 1. RAYLEIGH-RITZ RESULTS AT CONTACT POINT 


VS | Current Model Analvtical Solution | 
MPa) 285. 
esa 






a2 





Fieure 10. 


1.0 L.2 


NORMALIZED STRESSES 
0.6 


0.2 


0.0 


Comparison 


0.5 


1.0 1.5 2.0 
NORMALIZED DEPTH BELOW CONTACT SURFACE 


of o, in a roller contact problem 


oS 





1.0 


NORMALIZED STRESSES 
0.8 0.8 


0.4 





Figure 11. Comparison of o, in a roller contact problem 
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Figure 12. Stress contour of o, from Rayleigh- Ritz method 


2. Finite Element Method Results 
[In Chapter [II, an approximation method was developed to solve the contact 
problem using the finite element method. A method of calculating a svstem’s strain en- 
ergy and the total potential energy was investigated. In addition, the use of static 
condensation to improve optimization efficiency was described. As discussed, the nu- 
merical minimization techniques Were again utilized to determine the equilibrium posi- 
tion. A means of employing these techniques to evaluate the contact phenomenon was 
introduced. In this section, a simple contact problem will be investigated to validate the 
algorithms used to calculate strain energy and those used to implement static 
condensation. With confidence in these algorithms, the more complex roller-foundation 
problem will be examined. 
a. Iwo Thin Plates in Contact 
The procedure developed was first validated on a simple contact problem 


the solution of which was known. In this example, two thin bodies in plane stress were 
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Figure 13. Stress contour of o, from Rayleigh-Ritz method 


studied. As shown in Figure 16, one body, restrained on one edge and subjected to a 
horizontal load, comes in contact with a second body rigidly supported on three sides. 
The finite element model developed to solve this problem is composed of 14 linear tri- 
angular elements as shown in Figure 17. The objective of the fortran program developed 
to solve this problem was to calculate the total potential energy of the system using the 
finite element technique and the method of static condensation. With this accomplished, 
the equilibrium position can then be found via the augmented Lagrange multiplier 
method. 

The objective function for this problem is the total potential energy, 
Equation 3.14. Referring to Figure 16 and 17 the constraints imposed on the system are 


expressed as: 


ug <= 9.005 + u,, 
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Figure 14. Strain contour of ¢, from Rayleigh-Ritz method 


Uy < 0.005 + uy, 


where uw represents the horizontal deformation of node /. 


As shown in Figure 17, seventeen nodes were used to model the svstem. 
Each node has been assigned horizontal and vertical deformation vanables. Accord- 
ingly, the degree of freedom for this system 1s twice the number of nodes. As discussed 
in Chapter III, static condensation requires the identification of essential nodes and 
non-essential nodes. Essential nodes are those nodes where ADS design variables are 
assigned and boundary conditions are applied. Referring to Figure 17, node 3 is the 
point of load application and must be assigned a design variable. Nodes 8, 9, 11, and 
12 are assigned design variables in order to define the constraint equations described 
above. After eliminating all fixed nodes from consideration, nodes 2, 5, and 6 are the 


only candidates for condensation. As discussed in Chapter III, the global stiffness ma- 


a0 
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Figure 15. Strain Contour ¢, from Rayleigh-Ritz method 


trix 1s rearranged according to Equation 3.34. For this problem, the vector containing’ 


the non-essential nodal information, u,, is arranged as follows: 
ia 
{up} = {uy v2 Us Vs Ug V6} 


Where u, and v, represent the horizontal and vertical deformation of node /, respectively. 
Strain energy and total potential energy were calculated according to Equations 3.38 and’ 
3.14. The latter was minimized using the augmented Lagrange multiplier method de- 
scribed in Chapter II. 

The solution obtained from this simple problem were compared with the 
results obtained from Y.W. Kwon and J.E. Akin [ Ref. 15 ] and are shown in Table 2. 
The solutions were in agreement. It was concluded that a satisfactory procedure was in 


place to solve the more complex roller-foundation problem. 
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Figure 16. Two thin plates in contact 


Table 2. FINITE ELEMENT RESULTS: TWO PLATES IN CONTACT 


1 .Ox 10° 503x10-2 503x107? 
336x107 319x10— 


1.0x 10° 102x107? 134x10-? 
201x107 2235x1077 





b. Roller-Foundation Contact Problem 
A finite element gnd composed of 512 linear tnangular elements was con- 
structed to model the roller-foundation assembly shown in Figure 1. Because of the 
symmetry of the problem, one-half of the foundation was modeled. A refined mesh was 
constructed in the vicinity of the point of contact. The mesh 1s shown in Figure 18 


where the origin represents the point of contact. The domain dimensions are similar to 
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Figure 17. Finite element mesh for two plates in contact 
those chosen in the Rayleigh-Ritz method discussed in Part A of this chapter. The roller 
radius was chosen as 75 meters and the half-domain dimensions are 1.40 x 1.40 meters. 
As discussed in Part A, these dimensions represent an analysis of the region immediately 
adjacent to the contact zone and are a result of the local nature of the contact problem. 
Constraint equations were constructed according to the discussion of 


Chapter II Part B. Boundary conditions were imposed in the following manner: 


¢ Horizontal and vertical deformations were prohibited on the remote mesh bound- 
aries (1.e., u( 1.40, y) = (1.40, y) =0 and u(x,1.40) = v(x,1.40) = 0). 


e Horizontal deformation was prohibited on the axis of symmetry (i.e., u(0, y) = 0). 
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Figure 18. Finite element mesh for roller-foundation problem 
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Since the system has 289 nodes, the resulting 578 degrees of freedom ne-| 


cessitated the utilization of static condensation. Referring to the contact surface, the 







deformation variables that correspond to nodes on this surface are required for con- 
straint equations. The nodes that comprise the other three borders of the domain are’ 
all subject to boundarv conditions. Consequently, the interior nodes of the mesh are the 
nodes that are candidates for stauc condensation. In this model, all interior nodes were. 
condensed. The original svstem was reduced from 578 to 128 degrees of freedom. After: 
application of the boundary conditions, there were 46 possible deformations, a suffi-: 
ciently small number of design variables for the optimization algorithm. One additional 
design variable was used to represent the distance of travel by the roller. This value is: 
needed to calculate the work done by the roller on the bearing foundation. 

As before, Equation 3.38 was used to calculate the svstem’s strain energy. 
Following calculation and the subsequent minimization of the total potential energy, a 
post-processing procedure was followed to determine the contact stresses. The output. 
of the optinuzation routine represents the nonzero components of the {u,} vector. In) 
order to calculate stresses throughout the body, the remaining deformations contained! 
within the condensed vector {u,} must be determined. This vector is calculated directly 
using Equation 3.35. With deformations known throughout the domain, strains can be 
determined by applying Equation 3.27 to each element. The subsequent application of 
Equation 3.18 enables determination of stress for each element. 


To illustrate the capability of this method, an isotropic material with the 


following properties was chosen: 


E = 240 GPa 
v= 0.3 
G = 9253°GRa 


Load = 90.0 MPa 


Figure 19 shows the deformation resulting from the loading. For clarity, the defor- 
mations have been magnified 100 times their original values. Comparisons of the stress 
distributions with the analytical solution along the axis of symmetry are shown in Fig- 
ures 20 and 21. These figures are similar to Figures 10 and 11 and represent normalized 
versions of contact stresses. As shown in these figures, this method is a good approxi- 


mation of the stress distribution in the foundation of a loaded roller bearing. If the mesh 
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was more refined near the contact zone and the domain extended further, the agreement 
between the numerical and analytical solutions would be better. Figures 22 and 23 rep- 
resent normal stress contours of this problem. Figures 24 and 25 show normal strains. 
Table 3 shows a comparison of the results of this model and the analvtical solution at 


a selected clement in the region of contact. 


Table 3. COMPARISON OF STRESSES NEAR THE POINT OF CONTACT 


ution 
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B. APPLICATION 


The preceding section illustrated that the contact problem can be accurately simu- 


lated using the methods developed in Chapter III. It 1s the objective of this section to 
show how this approach can be applied to a contact problem in a composite plate sub- 
jected to low-velocity impact. 

A multi-ply laminate model has been constructed to investigate the response of 
composite materials to low velocity impact. It has been found that composite bodies 
subject to impact damage commonly fail due to delamination. Sandwich composites are 
currently being considered for use as turbine blades. I[t would be beneficial to acquire 
an understanding of the behavior of sandwich materials to impact damage. 

[In order to accomplish this task, a clamped composite beam similar to the one de- 
picted in Figure 26 has been modeled. The beam length 1s 25 cm. Beam thickness ts 2.5 
cm. Because of symmetry, half the beam was modeled with 256 bilinear elements. As 
Spoewn in figure 27, the mesh is refined near the point of contact, the origin of the mesh. 

The major assumptions of this model are that the beam 1s in a condition of plane 
strain and that the dynamic effect of the impact can be neglected. Reference {0 ap- 
proximated the loading resulting from low velocity impact as a sinusoid. In this study, 
the peak load will be examined to study the maximum normal and shear stresses. Ac- 
cordingly, the stress distributions obtained from this study will represent a “snap-shot’ 
in time of the response of the body at the instant of maximum loading. 

The finite element program developed to solve this problem 1s sufficiently flexible to 


alter the material stiffness matrix [D] shown in Equation 3.18 during construction of the 
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Figure 19. Deformed finite element solution: deformation magnified 100 times 
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Figure 20. Analytical solution vs. approximate o, from finite element results 
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Figure 21. Analytical solution vs. approximate o, from finite element results 
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Figure 22. 
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Stress contour o, from finite element results (increment 0.01 GPa) 
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Figure 23. Stress contour o, from finite element results (increment 0.02 GPa) 
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Figure 24. Strain contour <, from finite element results (increment 0.000035) 
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Figure 25. Strain contour ¢, from finite element results (increment 0.00001) 


Figure 26. 


Clamped composite beam 
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Figure 27. Finite element mesh for clamped beam model 


finite element global stiffness matrix. Therefore by defining the lavup for the laminate, 
the lamina stiffness matrices can be varied from element to element to accurately model 
the behavior of the body. This enables a variety of laminate lavups and lamina thick- 
nesses to be studied. 

The sandwich materials used in this study are composed of an isotropic interior 
material and orthotropic exterior laminae. Since a condition of plane strain was as- 
sumed, the isotropic material stiffness matrix is given by Equation 3.19a. For the 
orthotropic exterior laminae, the material stiffness matrix is given by Equation 4.1] 


(Ref. 16 ]. 
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ioey any ee 
E,.= Young's Modulus in / direction 


v= Poisson's rauo for lateral contraction in /* direction resulting from loading in 


the ¢ direction 


The derivation of the total potential energy calculation in Chapter III, Part B was 
done using linear triangular elements. Since bilinear elements were used in this model. 
calculation of the element stiffness matrix was more computationally intensive. Indi- 
vidual entries of the element stiffness matrix were obtained from Reference 17. Other- 
wise, the calculation of total potential energy was identical to the procedure outlined in 
Chapter ITT. 

There were two groups of boundary conditions applied to the problem. Along the 
clamped edge, deformation was prohibited. In addition, horizontal deformation was 
prohibited along the axis of symmetry. As before, the application of static condensation 
requires the idenufication of essential nodes the information of which 1s contained within 
the {u,} vector. In addition to the essential nodes associated with the above boundary 
conditions, the nodes along the contact surface are needed for the constraint equations. 
All other nodes were condensed out. 

To illustrate this application of problem solving, sandwich material with the fol- 
lowing properties Were studied: 


Exterior Laminae 


Ey, = 170 GPa 
EB, =alleSsGRa 
Go =e). Gra 
Vse—aO oS 


Isotropic Core 


iE = 225Gra 
G = 08S9GPa 
v = 0.35 





The beam was loaded by contact with a 10 cm radius ball. The transmitted force) 
was 250 N. Some unexpected trends were observed in the equilibrium position deter- 
muned by the optimizer. By examining the deformations along the axis of symmetry, a 
gradually decreasing trend in deformation moving away from the point of contact was: 
interrupted in lamina of significantly decreased stiffness. It is believed that this difficulty 
resulted from an inability of the optimization routine to approximate deformations 
through regions containing verv different orders of strain energy. In spite of these dif- 
ficulties, some critical information was obtained from this program. -As discussed in the: 
beginning of this study, one of the greatest difficulties of the contact problem is the de- 
termination of the size of the contact zone. Fortunately, the size of the contact zone cam 
be readily determined by examining the output from the constraint equations. By com- 
paring the ball radius (r) and the distance between the ball center and the node (r’), it 
can be determined if a node 1s in contact. This condition is illustrated in Figure 28. The 


distance r’ to the 1,, node 1s given by the equation: 
’ / = 2 2 
r=. (r-—o+) +x; 


If r’ is greater than r, the node ts not in contact with the body. 

With the extent of the contact zone known, the solution to this problem was ob- 
tained bv applving the contact boundary condition to a direct finite element program. 
Since the validity of the optimization program was in question, this method was applied 
using a 0-90-0 lavup similar to one used in a study conducted by Sun and Rechak [Ref. 
10]. The solution obtained from the current approach was verv close to the other re- 
sults. 

With confidence in the procedure, it was desired to examine the behavior of this 
model to various layups. The objective was to illustrate how this solution technique can 
be used for meaningful research. The study conducted by Sun and Rechak analyzed 
methods of reducing the likelihood of composite failure due to delamination. Of par- 
ticular concern is the magnitude of shear stress distribution and tensile stress in the y- 


direction, the two predominant causes of delamination failure. Using materials with the 
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Figure 30. Stress contour: r,, for clamped beam model 


thickness of the exterior laver increases, the maximum value of this stress increases. 
However, this stress is alwavs compressive at the interface with the core. By comparing 
the magnitudes of the shear and transverse normal stresses at the interface, it would 
appear that if delamination was to occur at this interface, it is more likely to be caused 
bv high shear stresses. 

Comparisons of the bending stresses at cross section A are shown in Figures 39, 40, 
and 41. These figures show that as the thickness of the exterior layer increases, a non- 
linear stress distribution intensifies in the layer closest to the contact surface. This trend 
would indicate that beam theory 1s unsuitable for estimating bending stress through this 
lamina. This nonlinear behavior is local to the contact zone. Figure 42 shows the 
counterpart for Figure 4] at cross section C. The stress distribution in the contact laver 
IS approximately linear. Another significant observation can be made by examining the 
bending stress graphs. As the thickness of the exterior layer opposite to the contact 


laver increases, the stress distribution within this layer transforms from purely tensile 
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Figure 31. Stress contour: o, for clamped beam model 


behavior to compressive-tensile behavior. This is significant because a bending crack 
initiated by tensile stresses tends to propagate to the core interface and cause delami- 
nation. The presence of compressive stress within this laver will tend to slow the growth 
of this crack toward the core. 

Thus far symmetric layups have been studied. To analyze how beams with asvm- 
metric lavups respond to contact loading, two beams with the following designations 
were studied: 0(3)-[SO(6)-0(7) and 0(8)-I[SO(6)-0(2). The first designated layer is the 
lamina closest to the contact surface. Figures 43 and 44 show the shear stress distrib- 
utions for these two layups. As before, the maximum shear stress 1s relatively unaffected 
by the different layups. As was the case for the symmetric beams, a thicker exterior laver 
close to the contact zone results in a more gradual transition of shear stress into the 
core. The result is a lower shear stress at the interface for the 0(8)-I[SO(6)-0(2) case. 
As seen in Figures 45 and 46, the trends for the transverse normal stress, o,, Were the 


same as those found in the symmetric beams. Deflection was less for the 





Figure 28. Determination of nodal contact 


properties outlined above, a sandwich composite beam with outer fibers aligned to 0 
degrees was first studied. Since the finite element model 1s composed of 16 lavers, the 
beams studied will be descnbed with the number of finite element lavers in parenthesis 
following the laver description. For example, a 0(3)-1SO(10)-0(3) beam 1s composed of 
10 isotropic core layers within 3 lavers of maternal with the fibers onented at O degrees 
on the top and bottom of the beam. 

Three symmetric layups of varying core thickness Were initially studied. The beams 
have the following designations: 0(3)-1SO(10)-0(3),  0(4)-TSO($)-0(4), 9 and 
O(5)-ISO(6)-O(5). Figure 29 shows the deformed 0(3)-1SO(10)-0(3) beam with defor- 
mations magnified 100 times. Stress contours for this beam are shown in Figures 30, 
31, and 32. Figure 30 shows the shear stress contour for the loaded condition. This 
stress 1s of particular concern since delamination, a common failure mode for compos- 
ites, 1S Commonly initiated by high shear stresses or tensile transverse normal stresses. 
As the figure shows, a very high stress gradient 1s present near the contact Zone. As the 
distance along the beam increases away from the contact zone, the magnitude of the 
gradient decreases until the cross sectional shear stress distnbution becomes parabolic. 
The transverse normal stress is also concentrated around the contact zone. 

Figures 33, 34, and 35 show the cross sectional shear stress distributions for the 


three symmetric beams described above. Three separate cross sections are shown on 
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Figure 29. Loaded beam: deformation magnified 100 times 


each graph. the locations of which are ‘ndicated in Figure 27. The vertical dashed lines 
on each graph identify the lamina interfaces. These figures show that all significant ac- 
tivity is confined to the lamina closest to the contact Zone. It is also evident that the 
maximum shear stress 1s relatively insensitive to the core thickness. Flowever, as the 
thickness of the layer closest to the contact surface increases, the shear stress transiuon 
is more gradual at the interface with the core. The result is a reduction in the shear 
stress at the interface for layups with thicker exterior lamina. It is also noteworthy that 
the maximum shear stress occurs at cross section B vice cross section A. 

Figures 36, 37, and 38 are graphs for the transverse normal stress. 6, , for the same 
three symmetric layups. These graphs show an increase in a, as the thickness of the 
exterior layers increases. The increased thickness of these layers produces a strongel 
beam, hence beam deflection and contact zone size are reduced. Accordingly, contact 
stresses increase. At Cross section C, some tensile transverse stress is evident. As pre: 


viously stated, tensile transverse stress is a potential source of delamination. As th 
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STRESS-Y: ZOOM OF CONTACT ZONE 
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Figure 32. Stress contour: o, for clamped beam model (In region of contact) 


0(8)-1SO(6)-0(2) case. The resulting smaller contact zone lead to higher contact stresses. 
mvs scen in the svmiumetric cases, an increase in the thickness of the laver closest to the 
contact zone resulted in an increase in the magnitude of the maximum tensile transverse 
stress, seen at cross section C. However, the stress at the laminate interface was alwavs 
compressive. 

With regard to bending stresses for these lavups, Figures 47 and 48 clearly show the 
nonlinear behavior as the contact layer thickness increases. In addition, the thickness 
of the exterior layer opposite to the contact surface shows similar results as the sym- 
metric cases. As the thickness of this layer decreases, the beam is more susceptible to a 


bending crack that propagates into the interface. 


Figure 33. 
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Figure 34. 
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Figure 35. Stress distribution: t,, for 0(5)-ISO(6)-0(5) laminate 
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Figure 36. 
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Figure 37. Stress distribution: o, for 0(4)-ISO(8)-0(4) laminate 


Figure 38. 
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Stress distribution: o, for 0(5)-[SO(6)-0(5) laminate 
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Figure 39, 
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Figure 40. Stress distribution: o, for 0(4)-ISO(8)-0(4) laminate 
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Figure 41. Stress distribution: o, for 0(5)-I[SO(6)-0(5) laminate 
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Figure 42. Stress distribution: o, for 0(5)-[SO(6)-0(5) laminate at cross section C 
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Figure 43. Stress distribution: t,, for 0(3)-[SO(6)-0(7) laminate 
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Figure 44. Stress distribution: t,, for 0(8)-1SO(6)-0(2) laminate 
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Figure 45. 
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Figure 46. 
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Figure 47. 
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Figure 48. 
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Vv. CONCLUSIONS AND RECOMMENDATIONS 


This study has developed two methods for approximating contact stresses using the 


augmented Lagrange multiplier method. As illustrated in Part A of Chapter IV, these 


methods accurately approximate the stresses that result from contact between a cylinder 


and plane surface. This study has also illustrated how this approach can be applied to 


understand the behavior of an actual contact problem by examining the response of a 


composite plate to low velocity impact. 


A. 


RAYLEIGH-RITZ APPROACH 


In the process of developing these methods, a number of comments can be made 


regarding the application of the Ravleigh-Ritz method to solving contact stress prob- 


lems. 


l. 


The selection of the trial function is an extremely challenging process. If it is de- 
sired to determine the deformation in a contact problem, the proper stress field 
must be first satisfied. Because of this, the selection of possible trial functions 1s 
limuted. For example, when selecting a trial function for the vertical deformation 
of the foundation of Figure |, a suitable selection is given by the following 
equation: 





at — sie cos( = ) 


This equation exhibits the favorable characteristics of maximum deformation at the 
contact surface and diminishing deformation as the distance from the contact sur- 
face increases. If contact stresses are to be modeled, this trial function ts inappro- 
priate. Calculation of ¢, is as follows: 


Ov oe TY 
ey Sage ele * ( 2] 


This function exhibits zero strain at the point of contact increasing to maximum 
Strain at the lower boundary. 











Since the selection of trial functions is difficult, the task is further impeded by 
complicated geometries. Furthermore, selection of a trial function necessitates that 


some knowledge of the deformation field exists. Without a sensible selection of 


trial functions, an accurate approximation is unlikely. 


This method assumes the trial function in the form of an infinite series. Solution 
accuracy theoretically should improve with an increased number of terms. How- 
ever, precautions must be taken to ensure the solution is numerically stable as the 
number of terms increases. Since the strain energy calculations require integration, 
there are choices of trial functions that will increase without bound as the number 
of terms is increased. This problem can be controlled by normalizing dimensions 
or limiting the choice of trial functions. 
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4. An increase in accuracy was observed as the number of constraints was increased 


B. 


and the distance between consecutive constraints was decreased. It is believed that 
the improved accuracy results from a better definition of the contact surface. 


FINITE ELEMENT APPROACH 
The results in Chapter III illustrated that this approach of applying the finite ele- 


ment method to contact stress analysis is effective. A number of comments can be made 


regarding this approach to problem solving. 


le 


oe) 


C. 


‘As illustrated in the results, this method accurately approximated the isotropic 
roller bearing problem. However, some difficulties were encountered during the 
modeling of the multi-ply composite. [In this model, smooth trends of decreasing 
deformations were often interrupted by spurious deformations or groups of defor- 
mations. These interruptions occurred within layers of significantly reduced 
suffness. It is believed that the optinuzation routine had difficulty approximating 
the deformations through these lavers because of their very small contribution to 
Strain energy. -\s stated in the results, the contact boundary conditions were ob- 
tained from the optimization program and applied to a direct finite element pro- 
gram to solve the problem. The above difficulty is recognized as a limitation of this 
approach. 


This approach is much more flexible for complicated geometries than the 
Ravleigh-Ritz approach. In addition, detailed knowledge of the deformation field 
is not needed as required bv the Rayleigh-Ritz approach. 


The application of static condensation 1s crucial to the successful implementation 
of this method. Every effort should be made to reduce the number of design vari- 
ables to improve optimization efficiency. 


COMMENTS ON OPTIMIZATION 


A number of observations were made regarding the general use of the Automated 


Design Synthesis System and the specific usage of the augmented Lagrange multiplier 


method. 


|. 
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The global optimum was more likely to be determined when the objective function 
was normalized. 


For both the Ravleigh-Ritz approach and the finite element approach, the initial 
choice of design variables had a significant affect on the possibility of obtaining the 
global optimum. For the Ravleigh-Ritz method, initial selections of design vari- 
ables can result in largelv dissimilar values of strain energy and external work, the 
two components of the objective function. [t was determined that convergence was 
more likely when optimization commenced with these two terms on the same order 
of magnitude. With regard to the finite element approach, sensible choices of the 
initial design variable vector was necessary for convergence to the global optimum. 
This was accomplished by intuitive selection of design variables to model the likely 
deformation. 


Solution accuracy can be improved by scaling constraint equations. It has been 
stated that in some circumstances, some constraints change more rapidlv than 


va) 
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others and can influence the solution excessively while others have little influence 
[ Ret. 6: pelse) 


With regard to the usage of the augmented Lagrange multiplier method, it was 
frequently necessary to ‘tune’ the optimization algorithm to a specific problem. 
This was done by varying the initial penalty term p and the initial Lagrange multi- 
plier term 4. As stated by Vanderplaats, commencing optimization with a small 
value of p should theoretically suffice for most problems [Ref. 6: pp. 137-138]. 
However, it was frequently necessary to select an initial value for p due to conver- 
gence to unrealistic solutions. Similarly , an initial choice of the Lagrange multi- 
pher term can effect the solution. Commencement with a small value is again 
recommended. [Ref. 6: p. I41]. This need to ‘tune’ the problem is a significant 
drawback to using this optimization method. The ideal way to overcome this lim- 
itation 1s to first tune the optimization routine using a known solution. With this 
accomplished, this approach can be used for meaningful data collection. 


SANDWICH COMPOSITE MATERIAL STUDY 


The behavior of sandwich composite materials to low velocity impact loading was 


successfully investigated by the application of the finite element approach. A number 


of observations can be made from examining the results. 
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E. 


The maximum shear stress 1s relatively insensitive to laver thicknesses. However, 
as the thickness of the contact laver increases, a reduction of the interface shear 
stress 1s observed. 


Tensile transverse normal stresses exist at some cross sections away from the con- 
tact zone. However, this stress 1s alwavs compressive at the interface. Compressive 
transverse stresses increase in beams with smaller cores due to reduced deflection 
dnd contact zome size. 


As the thickness of the layer closest to the contact zone increases, a nonlinear dis- 
tribution of bending stress within this layer intensifies. This phenomenon is local- 
ized to the region of contact. 


As the thickness of the layer opposite to the contact zone increases, bending crack 
propagation toward the core 1s less likely due to increased compressive bending 
stresses within the layer. 


RECOMMENDATIONS FOR FURTHER STUDY 


The methods developed in this study offer a basis from which additional research 


can grow. A reasonable direction is the relaxation of some of the assumptions made in 


Chapter If Part B. For example, relaxation of the rigid roller assumption and the 


frictionless surface assumption would provide challenging research. Models with com- 


plex geometry could be created. For example, a model of a pin loaded bolt connection 


could be created with ngid or non-rigid pins. Implementation of these changes would 


provide a versatile and highly applicable model for contact stress analysis. 
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